Crystallization and melting of bacteria colonies and Brownian Bugs 
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Motivated by the existence of remarlcably ordered cluster arrays of bacteria colonies growing in 
Petri dishes and related problems, we study the spontaneous emergence of clustering and patterns 
in a simple nonequilibrium system: the individual-based interacting Brownian bug model. We map 
this discrete model into a continuous Langevin equation which is the starting point for our exten- 
sive numerical analyses. For the two-dimensional case we report on the spontaneous generation of 
localized clusters of activity as well as a melting/freezing transition from a disordered or isotropic 
phase to an ordered one characterized by hexagonal patterns. We study in detail the analogies 
and differences with the well-established Kosterlitz-Thouless-Halperin-Nelson- Young theory of equi- 
librium melting, as well as with another competing theory. For that, we study translational and 
orientational correlations and perform a careful defect analysis. We find a non standard one-stage, 
defect-mediated, transition whose nature is only partially elucidated. 

PACS numbers: 02.50.Ey, 05.10.Gg, 05.40.-a, 87.18.Ed 



INTRODUCTION 

Spontaneous emergence of clustering is a widespread 
phenomenon in population biology, ecology, material sci- 
ence, and other fields [l|, Either passive or active 
"particles" (trees, bacteria, plankton, etc.) bunch to- 
gether forming dense localized clusters embedded in an, 
otherwise, almost empty surrounding space. The patchy 
distribution of plankton in the ocean surface Q , the spa- 
tial distribution of vegetation in semi-arid regions Q , or 
the fascinating patterns generated by bacteria grown in 
Petri dishes [5|,|6| are some examples. Some simple mech- 
anisms leading to clustering have been described: 

(i) Non-interacting inertial particles moving in a fluc- 
tuating eiiyironment (as a turbulent flow) may become 
clustered [7] . This path coalescence mechanism has been 
claimed to contribute to plankton patchiness. 

(ii) Branching and annihilating Brownian particles 
(also called super- Brownian processes) tend to bunch to- 
gether. This type of clustering has its roots in the fact 
that offsprings are created within a local neighborhood of 
their parents and die anywhere, giving rise to an overall 
tendency to form localized colonies [8|]. 

(iii) In some cases, clustering stems from the inter- 
action with a second "field" influencing the dynamics 
(water in models for vegetation patchiness 01, nutrients 
or chemicals in models for chemotactic bacteria cluster- 
ing 0, etc.) through some type of feedback mechanism. 
Reaction-diffusion equations have proven to be a conve- 
nient tool to model clustering at a mesoscopic level in 
this case 



Contrarily to the cases above, there are clustering mod- 
els in which individuals interact with each other in a di- 
rect way. For example, their effective interaction can be 



such that the reproduction rate (or some other individual 
dynamical property) is diminished by a factor depend- 
ing on the relative abundance of neighboring individuals. 
As a well documented example of this, let us mention 
the empirically observed J anzen- Council effect in ecol- 
ogy: owing to various circumstances (presence of special- 
ized pests or predators, competence for resources, etc.), 
the effective reproduction rate of an individual decreases 
with the number of conspecific specimens surrounding 
it Jl^. Similar effects appear also in bacteria colonies, 
social phenomena, and in systems exhibiting collective 
motion 

In many of these systems, clusters are distributed in 
space in a disordered way but, very remarkably, in some 
other cases they self- organize forming rather regular pat- 
terns. As a simple example to bear in mind, let us 
focus on the strikingly ordered patterns self-generated 
by Escherichia Coli, Salmonella typhimorium, and other 
bacteria grown in Petri dishes. When grown in a sub- 
strate of nutrients under adequate conditions these bac- 
teria colonies structure themselves into clusters which, in 
their turn, self-organize in spirals, squares, or crystal-like 
hexagonal arrays (as a beautiful illustration, see Fig. 1 in 
or page 529 of [i). 

How ordered can such two-dimensional patterns be? 
This question resembles very much the problem of crys- 
tal ordering in equilibrium solids and the existence of a 
melting (freezing) solid/liquid (liquid/solid) transition. 
For two-dimensional systems in thermodynamic equi- 
librium the Mermin-Wagner-Hohcnberg (MWH) theo- 
rem rules out continuous symmetries to be spon- 
taneously broken in the presence of fluctuations. This 
covers the case of translational invariance and, there- 
fore, two dimensional crystals/solids cannot exhibit true 
long-range translational order (destroyed by low-energy 
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Goldstone modes). Indeed, the most popular equilib- 
rium meltin g th eor y (to be described in some detail be- 
low 15, [igI. IITI. OI) predicts the melting transition to 
occur in two-stages, and to include an intermediate hex- 
atic phase in between a quasi long-range ordered solid 
phase and the disordered (or isotrorac or liquid) phase. 
Alternatively, a competing theory [l9| predicts a unique 
discontinuous transition between such two phases. 

Note, notwithstanding, that all the clustering prob- 
lems described above occur away from thermodynamic 
equilibrium. Actually, many of them exhibit absorbing 
states, a quintessential example of irreversible nonequi- 
librium dynamics (fluctuations can lead from activity to 
the absorbing state, but not the other way around [20|). 

Does the MWH theorem applies to nonequilibrium 
problems as self-organized bacteria colonies? Do they 
exhibit a melting/freezing transition similar to that oc- 
curring in equilibrium? Do they have an hexatic phase? 

Nonequilibrium problems violating the MWH theorem 
are well known. A notorious example is provided by 
the model of self-propelled particles proposed by Vic- 
sek et al. [isj which illustrates how communities (flocks, 
schools) of interacting individuals (birds, fishes) move co- 
herently, with true long-range order, adopting a collective 
non-vanishing velocity (which breaks the continuous ro- 
tational symmetry). Soon after the introduction of such 
a model it was proved in an elegant way that, indeed, true 
long-range order emerges owing to inherently nonequilib- 
rium effects [ilj (see also 22 1 for a related, though dif- 
ferent, nonequilibrium mechanism leading to true long- 
range two-dimensional ordering). Is the MWH violated 
in a similar way for the clustering problems described 
above? 

Aimed at answering the previously posed questions and 
to scrutinize the role of fluctuations in nonequilibrium 
clustered systems exhibiting self-organized patterns, we 
analyze here a simple model of individual interacting par- 
ticles 2^ 24 , 2^ exhibiting clustering as well as ordered 
patterns. 

Our main conclusions are as follows: while we find 
a solid phase with quasi long-range translational order 
and a disordered phase in analogy with equilibrium situ- 
ations, we do not find any evidence either of the hexatic 
phase predicted by the mo st p opular theory for equilib- 
rium melting [3, 0, 12, EM; '^O'^ S'^d the discon- 
tinuous transition predicted by competing theories [l9| . 
Instead, we report on a one-stage continuous transition 
analogous to what was reported in numerical investiga- 



ogy: emergence of clusters and ordering into hexagonal 
patterns. Afterwards, we report on its numerical analy- 
sis paying attention to the melting transition (including 
a careful analysis of topological defects) and compare our 
results with standard equilibrium melting theories. Fi- 
nally, we discuss our findings from a general perspective 
and present the conclusions. 



THE INTERACTING BROWNIAN BUG MODEL 
AND ITS LANGEVIN REPRESENTATION 

The individual-based model we study here, the "inter- 
acting Brownian bug" (IBB) model, was recently intro- 
duced by two of us [231, |2J|. It consists of branching- 
annihilating Brownian particles (bugs, bacteria) which 
interact with each other within a finite distance, R [l^. 
Particles move off-lattice in a d-dimensional [0, L]'^ inter- 
val with periodic boundary conditions, obeying a dynam- 
ics by which particles can: 

• diffuse (at rate 1) performing Gaussian random 
jumps of variance 2D, 

• dissappear spontaneously (at rate (3o), 

• branch, creating an offspring at their same spatial 
coordinates with a density-dependent rate A: 



X{j)^max{0,Xo-NRij)/Ns}, 



(1) 



where j is the particle label, Aq (reproduction rate 
in isolation) and (saturation number) are fixed 
parameters, and Nfifj) stands for the number of 
particles within a radius R from j. 

The control parameter is fi = Xq — (3o while the function 
maxO enforces the transition rates positivity. See [2^ 
where similar models in which the death or the diffusion 
rate are density-dependent are studied. 

The main phenomenology of the IBB is as follows 
[23I 24, [2^. For large values of ^ there is a station- 



tions of other equilibrium [26| and nonequilibrium 27 1 
melting problems. An explanation justifying our find- 
ings, supported by the analysis of topological defects and 
likely to be valid for other systems, is proposed. 

The paper is structured as follows. First we define 
the model, describe its basic ingredients and introduce 
a Langevin equation capturing all its relevant traits at a 
continuous level. We discuss briefly its main phenomenol- clusters start filling the available space they self-organize 



ary finite density of bugs (active phase) while for small 
values the system falls ineluctably into the vacuum (ab- 
sorbing phase). Separating these two regimes there is a 
critical point at some value /ic, belonging to the directed 
percolation (DP) universality class [25| characterizing in 
a robust way transitions into absorbing states [2^ . 

In the active phase, owing to the local-density depen- 
dent dynamical rules particles group together forming 
clusters (see Fig. [T] (a) and (b)) provided that the dif- 
fusion constant is small enough (for large values of D, 
homogeneous distributions are obtained). Such clusters 
have a well-defined typical diameter (which can be ana- 
lytically estimated [29|) and a characteristic number of 
particles within, which depend on the parameters R, Ns, 
and /i [2^, 24 1 . Well inside the active phase, when the 
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FIG. 1: Upper panels: Snapshots of two-dimensional IBB 
model in its stationary state (time t — 2 x IC Monte Carlo 
steps) in the: (a) disordered active phase, /i = 1.0, and (b) 
the ordered/patterned active one fi — 2.0; other parameter 
values are: Ns — 50, R = 0.1, D — 10~®. Lower panels: 
Snapshots of the two-dimensional Langevin equation, Eq.((2]) 
in its stationary state for n = 2.0, D — 0.25, R = 10, Ns — 
50 in the: (c) disordered active phase a — 2.5 and (d) the 
ordered/patterned active one, a — 1.5. Observe that the 2 
snapshots to the right (b and d) have a crystal-like hexagonal 
ordering, absent in the other two (a and c) . Also, the crystal- 
like packing is more evident in the continuous model for the 
chosen parameters, but no qualitative difference exist between 
the upper and the lower panels. 



in spatial structures with remarkable hexagonal order 
(see FigH] (b)). 

From the theoretical side, an important breakthrough 
is that the IBB model can be cast into a continuous 
stochastic equation 1231. Indeed, by applying standard 
Fock-space (Doi-Peliti jSOl) techniques, a Langevin equa- 
tion including the main relevant traits of the problem in 
a parsimonious way can be derived (see Appendix A in 
[23j and references therein). The Langevin equation for 
the local density [sll of bugs yo(x, i) reads (in the Ito 
representation) 



dpjy^, t) 
dt 



p(x, t) 



(2) 



x-y|<fl 



rfyp(y, t) + crVp(x, t)'q{x, t), 



where the noise amplitude cr is a function of the micro- 
scopic parameters and ?7(x,t) is a normalized Gaussian 



white noise. This includes only the leading terms in a 
density expansion; for instance, a higher order noise term 
appear in the mapping [2^ but it does not alter the re- 
sults reported in what follows in any significant way. 

Note that, leaving aside the non-local saturation term, 
this equation coincides with the Reggeon-field theory or 
Gribov process, describing at a coarse-grained level sys- 
tems with absorbing states in the directed percolation 
(DP) class (see (2^ and references therein). Let us un- 
derline the presence of a square-root multiplicative noise. 
Note also that the deterministic part of Eq.([2]), includ- 
ing a non-local saturation term, is identical to the equa- 
tion proposed by Fuentes et al. [llj. Eq.(I2) is therefore 
a simple (the simplest) stochastic generalization of such 
model. 

The main advantage of the continuous Langevin equa- 
tion above is that it allows for analytical studies and 
permits to scrutinize the effect of fluctuations (by just 
tuning the noise amplitude). It also constitutes a more 
general, elegant, and compact formulation of the prob- 
lem. For these reasons, we choose Eq.(I3]) as the starting 
point of our forthcoming analyses. 



INTEGRATION OF THE LANGEVIN EQUATION 
AND FIRST ANALYSES 



Analytical studies of the deterministic part of Eq.([21) 
(i.e. mean field analyses) have already been performed 
in [m, [2^ [l^] ■ They permit, for instance, to understand 
the wave-length instability mechanism leading to pattern 
formation, as well as many other aspects. In order to an- 
alyze the full stochastic Eq.Q and, in particular, its as- 
sociated melting/freezing transition, one needs to resort 
to computational studies. 

However, integrating numerically Langevin equations 
with square-root noise is a highly non-trivial task; ow- 
ing to the fact that for small density values the square- 
root term (with multiplies a random number) becomes 
larger in amplitude than the deterministic terms, stan- 
dard integration schemes (Euler, Runge-Kutta, etc.) lead 
ineluctably to unphysical negative values of the density 
field [13], and this pathology is not easily cured in any 
nai've way. Luckily, a very efficient scheme, specifically 
devised to overcome such a difficulty, has been recently 
proposed js^]. 

The method is a split-step algorithm in which the sys- 
tem is discretized in space and two evolution operations 
are performed at each discrete time-step: (i) first, the 
noise term is treated in an exact way, by sampling the 
conditional probability distribution coming out of the 
(exactly solvable) associated Fokker-Planck equation at 
each site. By sampling in an exact way such a distribu- 
tion an output is produced at each site, (ii) Afterwards, 
the remaining deterministic terms are integrated using 
any standard scheme taking as the input at each site the 
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output of the previous step at each site. More details and 
applications can be found in (33| . 

To implement the split-step scheme to integrate Eq.Q 
in two dimensions, we discretize the space, by introducing 
a lattice of linear size L — 256 or 512 [L « 1024 is al- 
ready at the limit of our present computational power). 
We fix the discrete time step to At = 0.25, R = 10, 
Ns = 50, D = 0.25 (which is small enough to have clus- 
tering) and use either ^ or cr as a control parameter. For 
all of the simulations reported here we initialize the sys- 
tem with a homogeneous initial density, p(x, t = Q) = pq, 
leave it relax towards its stationary state (reached typi- 
cally after 10^ Monte Carlo steps), and perform steady- 
state measurements (averaging over, at least, 10^ config- 
urations). 

First of all, we verify that Eq.(l2|) reproduces qualita- 
tively all the basic phenomenology of the microscopic IBB 
model: (i) an absorbing phase for small values of /i, (ii) 
an active disordered phase, encountered by increasing /i 
for a fixed a and (iii) an active ordered crystal-like phase 
which is reached by further increasing the value of /x or 
alternatively, fixing /i and reducing the noise amplitude, 
cr (see Fig.(5](c) and (d))). 

Separating (i) and (ii) there is a directed-percolation- 
like phase transition, while our focus here is on the tran- 
sition from the disordered active phase (ii) to the ordered 
self-organized one (in). In order to study the effect of the 
noise on such a transition, from now on, we fix ^ = 2.0, 
well into the active phase, and use the noise amplitude, 
tj, as a tuning parameter. Fig. [1] (plots (c) and (d)) 
shows two snapshots obtained for = 2.0 with a = 2.5 
and cr = 1.5 respectively. While in both cases clusters 
of localized activity (/o(x) ^ 0) exist, only in the second 
one clusters are self-organized into an ordered hexagonal 
array. 



Cluster analyses 

A preliminary step towards a systematic cluster anal- 
ysis is to have an efficient method to detect and label 
them. We have implemented an algorithm, based in the 
Hoshen-Kopelman one [sJl as follows. First, to avoid 
spurious clusters we apply a smoothening filter to the 
noisy field /5(x, t) in our simulations, which removes its 
short wavelength fluctuations: 



Pp(x,t) = 



1 



\y\=Rv 



p(x + y,t),dy 



(3) 



|y|=o 



where Rp is the cluster average radius (almost constant 
for the considered parameter range). Then we define a 
"smoothened" binary field, pp{x,t) taking a value (1) 
wherever Pp(x) < 9 (pp(x) > Q). The threshold Q is 
fixed to an optimal value O = 0.55, found by trial and er- 
ror, and we verify that an excellent cluster identification 
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FIG. 2: Main: cluster mobility per unit time (boxes and left 
axis) and mean cluster center of mass velocity (circles and 
right axis) as a function of a for L = 256. Inset: Number of 
clusters per surface unit versus cr, for L — 256. 



(as compared to visual inspection) is obtained. The out- 
put is not very sensible to the threshold choice although 
small variations can be observed. Once the binary dis- 
cretizcd field has been constructed, a standard Hoshen- 
Kopelman algorithm is straightforwardly employed. It 
generates a list of clusters, together with the list of spa- 
tial coordinates ascribed to (the center of mass of) each 
of them. 

Having described the main computational techniques, 
we now report on various observables characterizing col- 
lective properties of clusters. 

The average cluster mobility, m, is defined as the stan- 
dard deviation of the excursion of the cluster center of 
mass (during a fixed- length time interval), averaged over 
many different clusters in the steady state 



m 



^(^{x,-{x,)f) + ({y,-{y,)f), (4) 



where (xc, j/c) are the center of mass cluster coordinates 
and (■) stands for averages in the steady state. 

In Fig. [2] we plot m as a function of the noise ampli- 
tude cr; the mobility exhibits a sharp transition (sharper 
upon enlarging the system size) from the low-noise phase 
in which the clusters are almost frozen and localized in 
space to the high noise one in which they move more 
freely. The change of behavior occurs around a k, 2.39. 

Fig. [2] shows also the average velocity (in modulus) 
of the cluster center of mass [vx^T^y^) versus a. While 
for small values of a we observe a linear dependence be- 
tween {vc) and cr 35|, this linear dependence is broken 
above certain noise threshold (again around cr « 2.39), at 
which its derivative exhibits a discontinuity. Above the 
transition point the velocity increases non-linearly with 
cr. 
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We have also measured (Figl2] inset) the number of 
clusters per surface unit as a function of a. While in 
the disordered phase the density of clusters increases as 
the noise strength is reduced, it remains constant (at a 
value corresponding to the maximum capacity) once the 
threshold for an ordered structure is reached. 

The three described observables provide evidence of 
a melting/freezing transition. The change of behavior 
occurs in all cases at a unique point, somewhere around 
a w 2.39. A more detailed finite size scaling analysis 
would be required to pin down the critical point with 
more accuracy using these observables. 

The picture that emerges from these measurements is 
that clusters emerge at a mesoscopic scale out of the 
nonequilibrium microscopic rules and then, upon reduc- 
ing the noise-amplitude, they self-organize into frozen 
patterns with reduced mobility and velocity and with a 
more compact packing. In this sense, clusters become the 
equivalent of "particles" in standard liquid/solid transi- 
tions. 



Structure Function Analysis. 

In order to obtain an alternative, more direct, estima- 
tion of the location of the freezing/melting transition not 
relying on the (computationally costly) identification of 
clusters, we analyze a properly defined structure func- 
tion. As the overall density varies upon changing pa- 
rameters and system size, it is convenient to define a 
normalized version of the structure function as follows 



5(fc) = (|F(k)| 



(5) 



where F is the Fourier transform of the normalized den- 
sity 



F(k) = / /9„o™(x)e '•'■''dx. 



(6) 



k is the momentum, (■)^ stands for spherical averages 
over all two-dimensional vectors with module k = |k|, 
and the normalized density Pnorm is 



(x) = 



P(x) 



(7) 



Using this normalized density, and by virtue of the Par- 
seval's identity, /^^^ |F(k)|^ = j^^ |/Onorm(x)|^ = 1, it 
is guaranteed that S (fc) is normalized to unity for all 
parameter values and system-sizes. 

Fig [3] shows the structure- function for a size L — 512 
as a function of k, for three different noise amplitudes: 
a — 2.20, 2.40, and 2.60. The three curves exhibits a 
very pronounced first Bragg-peak at positions around 
koL « 43 (which corresponds to fco « 0.083(9) and 
therefore a separation between clusters Rq « 11.9(1). 
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FIG. 3: Main plot; Structure function, S{k), for three differ- 
ent values of the noise amplitude a and L = 512. The peak 
around feo ~ 42 becomes more pronounced as we go deeper 
into the ordered phase. Inset: S{ko) versus a for L — 256 and 
L = 512. 



Their respective heights decrease with noise strength: the 
larger the noise the less ordered the structure. Neverthe- 
less, perfect ordering (delta-peak at a given mode) will 
never be reached as cluster internal fluctuations enforce 
some dispersion in the structure function. 

Given that S{k) has been normalized, the height at 
fco, S{ko), constitutes a good measure of the degree of 
order. Actually, an abrupt transition is obtained at about 
a « 2.39 for the largest available size (see (FiglS])): above 
the transition S'(fco) decreases abruptly corresponding to 
the breakdown (melting) of the ordered patterns. Also, 
observe in the inset of FigO that the degree of order is 
enhanced upon enlarging the system size from L — 256 
to L = 512. We will come back to analyze this issue 
later. 



TWO-DIMENSIONAL MELTING THEORIES: 
A SHORT REVIEW 



The evidence accumulated so far, both from cluster 
analyses and structure function measurements, reveals 
the existence of a melting/freezing transition, somewhere 
around a = 2.39. The key question we face now is: 
does such a nonequilibrium transition exhibit the well- 
known universal features of standard equilibrium melt- 
ing/freezing in two-dimensional systems? 

Note, in addition, that here the number of clusters 
(particles) is not constant but fluctuates, and its average 
value varies as a function of the control parameter, a. 
Is this relevant for the properties at the melting/freezing 
transition? 

Before tackling these problems, for the sake of com- 
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pleteness, and for future reference, in this section we 
briefly summarize the main results of the celebrated stan- 
dard theory of two-dimensional melting: the Kosterlitz- 
Thouless-Halperin-Nelson- Young (KTHNY) theory [l5, 
[ItI . [lij . We also discuss briefly an alternative com- 
peting theory. 

The KTHNY is based on a statistical physics analy- 
sis of topological defects, i.e. particles with a number of 
nearest neighbors (assuming a Voronoi or Wigner-Seitz 
construction) other than 6. Dislocations perturb trans- 
lational order and disclinations hinder orientational or- 
der [13, ■ ^ detailed inspection of hexagonal ordering 
in the presence of fluctuations reveals that disclinations 
correspond to free monopoles, either five-fold or seven- 
fold, where five and seven refer to the number of near- 
est neighbors as measured in a Voronoi or Wigner-Seitz 
construction. Analogously, dislocations can be identi- 
fied with tight pairs (i.e. dipoles) of a five-fold and a 
seven- fold disclinations. The disordered (or isotropic or 
liquid) phase is characterized by the proliferation of de- 
fects: both monopoles and dipoles. 

The main prediction of the KTHNY theory is that, 
contrarily to what happens in higher dimensional sys- 
tems, where the melting occurs discontinuously at a 
unique transition point, in two-dimensional systems 
melting occurs in two stages. Translational and orien- 
tational order lose their stability at different Kosterlitz- 
Thouless-like jT^I critical points where dislocations and 
disclinations, respectively, unbind. 

The theory assumes that in the solid phase there are 
nor free dipoles nor monopoles, but only quadrupoles 
(low-energy excitations), that the number of disloca- 
tions/dipoles throughout the first (melting) transition 
is small and that they are generated progressively in a 
smooth way as the temperature is risen. This allows to 
treat the system at the melting transition as a weakly 
interacting gas of dipoles/dislocations. An analogous 
assumption is made for monopoles/disclinations at the 
second transition point where monopoles unbind from 
dipoles. 

The three phases put forward by the KTHNY theory 
are as follows 16, 17, [l3| (see Fig. [4] for a graphical 
illustration) : 

• Only at zero temperature the lattice ordering can 
be perfect while, for any non-vanishing tempera- 
ture, defects appear. Below a first critical point 
(denoted am here) defects are tight together in 
quadrupoles, hindering translational order. As a 
result, translational correlation functions decay al- 
gebraically in space (with continuously varying ex- 
ponents), as corresponds to quasi long-range or- 
der. Orientational correlations (see below for a pre- 
cise definition) decay at long distances to a non- 
vanishing constant value, corresponding to true 
long-range orientational order [s^ . 
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FIG. 4: Schematic presentation of the main predictions of the 
KTHNY theory for two dimensional melting. See the main 
text for more details. 



• Above a first critical point, am, 
dipoles/dislocations unbound from quadrupoles, 
destroy translational order (i.e. translational corre- 
lations decay exponentially fast). They also affect 
orientational correlations which decay algebraically 
with a continuously variant exponent tjq obeying 
J < '76 < ^ and a diverging correlation length, i.e. 
they exhibit quasi- long-range orientational order. 
This is the, so called, hexatic phase. 



• Above a second (isotropic) critical point, ai, 
monopoles/disclinations unbind from dipoles, hin- 
dering quasi long range orientational order. Both 
translational and orientational correlations decay 
exponentially. The associated correlation lengths 
diverge as stretched exponentials of the form 
aexp(6A~^/^) where A is the distance to ai. This 
is the isotropic (also called disordered or liquid) 
phase. 



This scenario has been verified in a number of numer- 
ical 17, [s^] and experimental [13, S studies, while it 



was not verified in others 

Despite of its success, the KTHNY is not the only 
plausible theory of two-dimensional melting. A com- 
peting one was proposed by Chui [l^, who argued that 
some systems should exhibit a unique first-order melting 
transition mediated by the appearance of "grain bound- 
aries" . In this picture, chains of defects appear limiting 
ordered grains and separate neighboring mismatching do- 
mains. The main difference between this theory and the 
KTHNY one is that here the transition appears owing to 
a collective excitation of defects. In some systems, it has 
been shown that the melting transition can change from 
KTHNY to first order upon changing some parameter 
(17l] . so both theories can be taken as complementary. 
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FIG. 5: Radial correlation function, g{r) as a function of r 
for different values of a and linear size L = 512. 




FIG. 6; Correlation length, ^ and exponent rj, obtained by fit- 
ting g{r) in Fig[S]to Eq.®. Note the abrupt increase of the 
correlation length around a = 2.39 from a small value above 
to a large saturating value (roughly equal to I//4 — 128, i.e. 
half of the maximum radius). The exponent rj changes con- 
tinuously, in remarkable agreement with the KTHNY predic- 
tions. 



ANALYSIS OF THE NONEQUILIBRIUM 
MELTING/FREEZING TRANSITION 

To determine the plausibility of the KTHNY scenario 
(or, instead, that of competing theories) for the nonequi- 
librium transition under study, we analyze in this section 
translational and orientational order in numerical simu- 
lations of Eq.llll). 

Translational Order 

Translational order can be studied by measuring the 
two-point radial correlation function: 

g{r = |r - r'l) cx (p(r)p(r')) - (p(r)) (p(r')) , (8) 

where (•) stands for averages in the stationary state taken 
over all pairs of particles at generic positions r and r' 
separated by a distance r = |r — r'|. 

In FigO we plot g{r) as a function of r for different 
values of a above and below the critical point. For 
all parameter values the wavy curves reveal the clus- 
tered nature of the density distribution. To determine 
their asymptotic trends, we analyze the envelope of 
such curves, and fit it to the behavior predicted by the 
KTHNY theory: 

g{r) oc r-"e-'-/« (9) 

where ij and ^ are fitting parameters. Note that, the 
power-law corresponds to the KTHNY prediction while 



the exponential allows to describe finite-size induced cut- 
offs. Figini shows the results of such a fit for different 
values of a. In all cases, the fit correlation coefficient is 
large than 0.92. Note first the abrupt jump of the trans- 
lational correlation length ^ at the critical value a w 2.39; 
in the ordered phases it converges to a saturation value 
controlled by system size (see figure caption), while in 
the disordered one it takes much smaller values. On the 
other hand, the exponent r/ changes continuously from 
a value nearby 0.31(1) at the critical point (in perfect 
agreement with the KTHNY prediction, which imposes 
1/4 < ri{a„i) < 1/3) to smaller values as a is decreased, 
confirming the generic algebraic decay of g{r) in the solid 
phase, with an exponential cut-off given by the system 
size. 

To check the internal consistency of our results, we 
also estimate r] by analyzing the previously obtained 
structure-function results. As S'(k) is trivially related 
to gir)) through a Fourier transform, and g(r) (which 
depends only on the modulus of r, r) can be modeled 
in the solid phase by g{r) ^ r^'' cos{kor) (for simplicity, 
we omit here the exponential, system size induced, cut- 
off e^*"/^) one can find after simple algebra, that for a 
two-dimensional system of size L 

S{k,L) ^ 271 [ r^-" cos{kor)Jo{kr)dr, (10) 
Jo 

where Jq is the zero-order Bessel function. Using the 
asymptotic behavior of Jo(fcr) it is not difficult to obtain 
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that 



^(fco,L)^L3/2-, 



(11) 



for the mode kg (which is the only one for which there 
are not destructive interferences). 

ActuaUy, in the inset of Fig. ^ we showed that while 
the curves for S{kQ,L) overlap for different system sizes 
above the critical point, bellow it, they growth with sys- 
tem size. Indeed they grow algebraically with a slightly 
tJ-dependent exponent which for all a values is in the 
interval [1.24,1.36]. Therefore, exploiting Eq. (fTT|) we 
derive a value for rj which is in the interval [0.14,0.26] 
in rather good agreement with the direct measurements 
above. Summing up, we have deduced, in two indepen- 
dent ways, results consistent with algebraically decaying 
correlations as predicted by the KTHNY scenario. 



Orientational Order 

To quantify the degree of orientational order of a given 
configuration in the steady state, we first identify the 
clusters using the algorithm described above. After that, 
a Voronoi tessellation of the system configuration is con- 
structed [3§|. For each configuration, the corresponding 
tessellation gives as output the list of clusters and the set 
of nearest neighbors of each. Finally, for each cluster j. 



we define [lEMIi 



^ k=l 



1 



(12) 



where the sum extends over the Nj neighbors of cluster j; 
9jk is the angle between the centers of mass of clusters j 
and k and an arbitrary fixed reference axis. The average 
of ■06 over different clusters, 



V'6 



1 ^= 

^E^6(j) 



(13) 



where Nc is the total number of clusters, is a global ori- 
entational order parameter. An associated susceptibility 

can be also defined, xe = Nc (^(V-'l) ^ (V'e}^) ■ 

Fig- El (main plot) shows ipe as a function of a for two 
different system sizes. For L = 512, V'6 changes abruptly 
in the interval from a large value for a below « 2.38 to rel- 
atively small ones above 2.40. Note that finite size effects 
operate in opposite directions below and above the transi- 
tion and the curves intersect at some point between these 
two regimes (resembling what happens for the Binder cu- 
mulant at continuous phase transitions). This provides 
a useful criterion to locate the critical point; using the 
two available sizes, the best estimation is cr = 2.39(2), 
and strongly suggests that the transition is continuous 
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FIG. 7: Main plot: Stationary orientational order parameter, 
ipe, as a function of a for two different system sizes, L — 256 
and 512. Inset: Orientational susceptibility, xe, as a function 
of (J for L = 256 and L = 512. 



in agreement with the KTHNY scenario. On the other 
hand, the susceptibility (inset of Fig. [7]) exhibits a sharp 
peak, which moves slightly to the right with increasing 
system size, being located at cr = 2.39(1) for L — 512. 

The probability distribution function (pdf) of 
provides additional information about the nature of the 
transition. In Fig. [5] we plot the pdf as a function of 
h/'6(j)l^ (to have only positive values) for various a and 
size L = 512. At small values of cr (as 2.00 or 2.35) 
the pdf is unimodal and peaked around a high value. 
Increasing the noise strength the peak becomes less pro- 
nounced. For noise values around the transition point, 
a e [2.40, 2.42] the pdfs are rather flat, a new peak ap- 
pears nearby zero, and the average value shift (in an ap- 
parently continuous way) from a high value in the ordered 
phase to a small one in the disordered one. 

To further check the predictions of the KTHNY theory 
it is necessary to determine the two-point orientational 
correlation function: 



96{r) = ('06(i)'06(fc)) 



(14) 



where the average is taken over all pairs j,k of clus- 
ters separated by a distance r (results shown in Fig. [5]). 
While for cr > 2.40 the envelope of the curves veers down 
in a double-logarithmic plot (i.e. decay exponentially) 
for values smaller than cr = 2.39 the envelopes converge 
to constant values in the long distance limit, signaling 
the emergence of true long-range orientational order in 
the crystal-like phase. We find no evidence of a hex- 
atic phase, characterized by generic power-law decaying 
orientational correlations. If it actually exists it should 
be very narrow and could not be detected within our 
present numerical resolution. Indeed, we have scanned a 
in steps of 0.005 units (not shown), and we always find 
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FIG. 8: Probability distribution functions of |i/;6(i)|'^ for dif- 
ferent values of a and linear size 512. 
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TABLE I: Correlation length and critical exponent rje at 
different values of a in the isotropic phase. Results obtained 
by fitting the (envelopes of the) curves for geir) to Eq.(|15[). 
The last column shows the correlation coefficient, Corr be- 
tween the numerical data and their fits. 
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FIG. 10: Orientational susceptibility, X6i and correlation 
length, ,^6, in the isotropic phase. The solid lines represent 
the best fit to the stretched exponential function, Ea. (|16|l . 



FIG. 9: Log-log plot of the (spherically averaged) orienta- 
tional correlation function, g6{r), for different values of a and 
linear size 512. 



that the envelopes either bend down below certain sepa- 
ratrix (i.e. the critical point, which here we estimate to 
be at cr = 2.395(5)) or flatten above such a value. 



Analysis of the isotropic phase 

As synthesized above, the KTHNY theory pre dicts 
that in the disordered or isotropic phase 3, 17, 3| 



(15) 



where tiq = 1/4 is a critical exponent and ^q, the ori- 
entational correlation length. The theory also predicts a 
divergence in as the critical point is approached 



e6(A) = a« 



(16) 



where A is the distance to the critical point and and 
6^ are constants. A similar, stretched exponential be- 
havior is also predicted for the susceptibility, X6(A) = 

Both ^6 a-iid can be measured by fitting the en- 
velopes of ge{r) for different values of a to Eq. (fT5|) . We 
have performed a two-parameter (^g and Tye) fit, and ob- 
tain the results summarized in Table HI Observe the very 
fast grow of upon approaching the transition point. 
The estimations of ryg are very close to the KTHNY value 
r]Q — 1/A (actually, they become indistinguishable by in- 
cluding error-bars). Not surprisingly, close to the critical 
point, the correlation coefflcient, Corr, is worse than far 
from it. 

In Fig. [To] we show the results obtained by fitting the 
simulation results to the predicted stretched exponen- 
tial behavior Eg. (fl6|) . The fit is excellent in both cases 
and allows for estimations of the critical point location: 
2.38(2) using and 2.37(2) using xe- 

These results provide a strong backing for a KTHNY 
scenario in the isotropic phase. 
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Summary of the numerical observations 

Summing up, the numerical analyses detailed above are 
consistent with the KTHNY theory in any respect, except 
for the fact that no evidence of a hexatic phase is found. 
Our results are compatible with a scenario very similar 
to the KTHNY one, but with the two transition points 
merging into a unique one, occurring somewhere between 
tJ = 2.38 and a = 2.39. No evidence is found of a hexatic 
phase. A scenario analogous to this, i.e. a one-stage 
continuous melting transition, has been found in other 
equilibrium and nonequilibrium systems 2^, 27 1. 



DEFECT ANALYSIS 

Inspired by the analyses in [27t . in this section 
we check the KTHNY assumptions rather than its 
predictions. As briefly explained above, the the- 
ory assumes that the number of dipoles/dislocations 
(monopoles/disclinations) throughout the first (second) 
transition is small and that they are generated pro- 
gressively in a smooth way as the temperature is 
risen. This allows to treat the first (second) transi- 
tion as a weakly interacting gas of dipoles/dislocations 
(monopoles / disclinations) . 

To check if this picture describes properly the behav- 
ior of our model, we scrutinize how defects appear and 
proliferate upon rising the noise amplitude. In Fig. [11] 
we plot defect maps for 4 different values of cr. The lines 
correspond to the Voronoi construction for a given con- 
figuration; five-fold, six-fold, and seven-fold clusters are 
marked with black circles, central dots, and red stars 
respectively (higher and lower order defects correspond 
to blank clusters). As the noise amplitude is increased 
the total number of defects grows. While below the crit- 
ical point only quadrupoles are significatively present, 
nearby the transition point isolated dipoles start un- 
binding. Even if unbound, they have a clear tendency 
to bunch together, and indeed, at slightly larger noise- 
amplitudes, in the isotropic phase, defects bunch together 
showing a tendency to form string-like structures. Deep 
into the isotropic phase defects proliferate, and extend 
through the system keeping, in any case, the propensity 
to form condensed string-like structures. 

In this respect, it is noteworthy that Fisher, et al. [40'| 
predicted, within the KTHNY theory, that monopoles 
can appear at correlated locations showing a strong ten- 
dency to arrange themselves in small-angle grain bound- 
aries. Contrarily, the condensates we find seem to be 
dominated by strings of dipoles. 

To quantify the phenomenology above, we measure the 
density of defects (quadrupoles, dipoles and monopoles) 
as a function of noise amplitude (see FigfT2|). The density 
of defects (of any type) grows from a small value below 
cr = 2.3 to a large value in the isotropic phase. The 



increase is very steep around cr = 2.38 for both dipole and 
monopole densities, suggesting the existence of a unique 
unbinding transition. 

To complement the previous plot, the inset of Fig. 
[T2l shows the density of isolated dipoles/dislocations and 
monopoles/disclinations versus cr, confirming the previ- 
ous observations, that the steepest increase occurs at 
cr « 2.39. Note that the total number of monopoles grows 
continuously upon increasing cr, while the presence of a 
relative maximum in the density of dipoles reflects the 
competition between two conflicting tendencies: the un- 
binding of dipoles from quadrupoles and the liberation 
of monopoles (and formation of condensed defects) from 
free dipoles. 

Finally, Fig. [13] shows the densities of isolated dipoles 
and of condensates as a function of the density of 
monopoles. This plot reveals a strong correlation be- 
tween monopoles and dipoles and, more importantly, it 
illustrates how they both vanish at roughly the same 
point, suggesting again (up to finite size effects) the ex- 
istence of a unique one-stage transition. Indeed, if the 
dipoles unbinding occurred before the monopoles one, the 
line joining the circles should intersect the vertical axis 
which is not the case. It is also at such a unique tran- 
sition point that condensed defects appear, being their 
density (roughly) linearly correlated with the density of 
monopoles. 

In summary, a detailed analysis of defects in Voronoi 
tessellations supports the interpretation that both 
dipoles and monopoles unbind at a unique critical point. 
At such a transition point condensed defects are gener- 
ated in an apparently continuous way. In a future pub- 
lication we will analyze in a more detail the correlations 
between isolated monopoles, dipoles and condensates to 
further clarify the analogies and differences with the stan- 
dard KTHNY scenario, and try to develop a theoretical 
framework for one-stage continuous melting transitions. 



DISCUSSION AND CONCLUSIONS 

Motivated by the observation of remarkably regular ar- 
rays of clusters formed by bacteria growing in Petri dishes 
and related problems, we have revisited the individual- 
based interacting Brownian bug (IBB) model in which 



birth rates are local-density dependent [23|, |2j|. Apart 
from an absorbing phase transition, this model shows a 
transition from a disordered active phase, in which par- 
ticles aggregate in localized clusters, to an ordered or 
crystallized active phase, in which clusters self-organize 
forming hexagonal arrays. For the sake of generality 
and aimed at facilitating analytical studies, rather than 
studying the discrete model itself we have analyzed its 
equivalent continuous Langevin representation. This 
stochastic equation, which involves a non-local satura- 
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FIG. 11: Defect maps at 4 different values of the noise amplitude, a and size L = 512. Lines define a Voronoi tessellation; 
five-fold, six-fold, and seven-fold clusters are marked with black circles, central dots, and red stars respectively. The (few) empty 
polygons correspond to 4-fold, 8-fold, and 9-fold defects, (i) Upper left panel (low noise intensity; a = 2.30); only localized 
quadrupoles appear, (ii) Upper right panel {a = 2.38 nearby the critical point); a proliferation of isolated dipoles/dislocations 
can be observed and a small number of isolated monopoles/disclinations can also be observed, (iii) Lower left panel (cr = 2.41 
slightly above the critical point); although isolated dipoles and monopoles exist, there is a clear tendency to form defect 
condensates with string-like geometry, (iv) Lower right panel (cr = 2.60 well above the transition point); defects proliferate and 
the tendency to form string-like structures is maintained. 



tion term as well as a square-root multiplicative noise, 
can be numerically integrated by discretizing it and em- 
ploying a recently introduced efficient integration scheme 
specifically designed to deal with square-root multiplica- 
tive noise [11]. 



The main results we obtain are: 

(i) First, we have shown explicitly that the continu- 
ous model (truncated to include only the leading terms) 
reproduces the phenomenology of the original discrete 
model, including an absorbing phase in which the sta- 
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FIG. 12: Main: Accumulated density of defects of the various 
possible types as a function of the noise amplitude, a. Inset: 
Relative isolated dislocation and disclination concentrations 
versus a. 
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FIG. 13: Dislocation density (left axis and circles) and total 
defect density (right axis and boxes) in logarithmic scale, as 
a function of the isolated disclination density {L — 512). 



tionary particle density is zero, a disordered active phase 
in which the density field is localized in clusters of activ- 
ity surrounded by empty regions and, finally an ordered 
phase with hexagonal patterns. 

(ii) We have studied the ordering transition by 
analyzing cluster properties (average velocity, mean- 
displacement, etc.) as well as by means of analyses of 
the structure-function. These studies reveal that a melt- 
ing/freezing transition indeed occurs at some value of the 
noise amplitude: for small noises, clusters are trapped 
into hexagonal configurations as a result of a collective 
effect, while for noise amplitudes above threshold they 



have much more mobility. 

(iii) To better understand and quantify the transition 
we have studied both translational and orientational cor- 
relation functions. In particular, we have verified that, 
in the solid phase, translational correlation functions ex- 
hibit generic power-law decay even if with system-size 
induced cut-offs. The value of the decay exponent at the 
melting transition being in excellent agreement with the 
KTHNY prediction. However, we do not found any ev- 
idence of an hexatic phase, contrarily to what predicted 
by the standard (Kosterlitz-Thouless-Halperin-Nelson- 
Young (KTHNY)) theory. This same conclusion is also 
borne out by analysis of the global orientational order- 
parameter. Up to the numerical resolution limit, this 
suggests the existence of a unique continuous transition 
point. 

(iv) Above such a transition point, i.e. in the disor- 
dered or isotropic phase, the correlation functions have 
been shown to decay as stretched exponentials, in excel- 
lent agreement with the KTHNY predictions. 

(v) We have performed a detailed analysis of de- 
fects and found that both free dipoles/dislocations 
and monopoles/disclinations seem to unbind from 
quadrupoles at a unique transition point, above which 
also defect condensates are formed. This is in contrast 
with the KTHNY scenario and leads to a continuous one- 
stage melting transition. 

(vi) The unique transition is continuous, so it cannot 
be explained by the main alternative theory to KTHNY 
[lit which predicts a first-order melting transition. 

Summing up, the phenomenology of this nonequi- 
librium model can be only partially described by the 
KTHNY theory. While the melting/freezing transition 
is indeed characterized by the smooth unbinding of de- 
fects this occurs through a unique continuous transi- 
tion. The reason for this seems to yield in the fact 
that once dipoles/dislocations unbind, the perturbation 
they generate around them is large enough as to un- 
bind also monopoles/disclinations: dipoles, monopoles 
and the string-like condensates they form are strongly 
correlated. 

Note that, the nonequilibrium microscopic dynamics 
of the interacting Brownian bug model (or its equivalent 
Langevin representation) is responsible for the generation 
of mesoscopic clusters. Once such clusters are generated, 
they interact in an effective way, and the physics at a 
macroscopic scale does not seem to differ in any essential 
way from other equilibrium problems [2^, \2l\ for which a 
similar one-stage continuous melting transition has been 
reported. We believe that this scenario is not specific of 
nonequilibrium systems but is determined by the way in 
which defects interact among themselves. This might not 
depend on the equilibrium versus nonequilibrium nature 
of the process but rather on other structural details in- 
fluencing the way defects interact among themselves. A 
more detailed analysis of defects and defect-correlation 
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will be investigated in a future work, where we will try 
to develop a theoretical framework for one-stage contin- 
uous melting transitions. 

Let us also emphasize that the models we have ana- 
lyzed are not the best choice to explore with high numer- 
ical resolution the possibility of one-stage melting from a 
general perspective. Effective models, with a dynamics at 
the level of clusters (as opposed to having a microscopic 
dynamics for particles) would be a much better option 
from the computational point of view. 

In summary, despite of interesting and not fully un- 
derstood differences, the striking patterns produced by 
the biologically inspired Langevin equation ([2]) resemble 
very much the melting/freezing solid/liquid transition in 
equilibrium systems. 

It is our hope that this paper will motivate further 
studies of (i) the effect of fluctuations on self-organized 
noncquilibrium patterns and (ii) the analogies and differ- 
ences between the defect-mediated type of melting tran- 
sition described here and standard equilibrium melting 
scenarios. 
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